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Abstract: We investigate the observational consequences of scalar instabilities in bimetric 
theory, under the assumption that the Vainshtein mechanism restores general relativity 
within a certain distance from gravitational sources. We argue that early time instabilities 
have a negligible impact on observed structures. Assuming that the instabilities affect 
sub-horizon density fluctuations, we constrain the redshift, Zi, below which instabilities 
are ruled out. For the “minimal” /3\ -model, observational constraints are close to the 
theoretical expectations of Z{ ~ 0.5, potentially allowing the model to be ruled in or out 
with a more detailed study, possibly including secondary cosmic microwave background 
constraints. 
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1 Introduction 

General relativity is the unique theory of a massless spin-2 field interacting with matter. 
While it has been tested to great precision within the solar system, its application on the 
level of galaxies and beyond points towards the presence of major unknown components of 
the universe—dark energy and dark matter—which have so far only been probed gravita¬ 
tionally. From a theoretical point of view, the observed value of the cosmological constant 
seems to require a high degree of fine-tuning in order to be reconciled with the quantum 
properties of the vacuum. Adding to this observational conundrum is the fact that gen¬ 
eral relativity is non-renormalizable as an effective field theory. These observational and 
theoretical problems have sparked a great deal of interest in possible extensions of general 
relativity in order to fully understand the nature of gravitational interactions. 

Following the classification of the particles in the Standard Model according to their 
mass and spin, a natural extension of general relativity is to introduce another spin-2 
field into the gravitational sector. This gives rise to both massive and massless spin- 

2 interactions. Such proposed interactions have a long history [1, 2], but finding the 
adequate theoretical formulation has proven difficult due to the presence of ghost modes 
in the spectrum [3]. Recently, a formulation, dubbed the Hassan-Rosen theory, with two 
interacting spin-2 fields was found that does not contain a ghost [4, 5]. This formulation 
contains the so-called de Rham-Gabadadze-Tolley theory [6, 7] in the limit when one of the 
fields becomes non-dynamical. That these theories were ghost free was proven in [5, 8-10]. 

The background cosmology of the Hassan-Rosen theory has been extensively studied 
[11-16]. It was shown, however, in Ref. [17] that certain instabilities could arise in the 
scalar sector on a perturbative level. The scalar sector was further studied in [18-30]. In 
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particular, Ref. [26] gave conditions for the existence of the scalar instabilities in terms of 
the parameters of the Hassan-Rosen theory. 

The instabilities are also present in the tensor sector, and have been studied in Refs. [31- 
34]. In all of these studies, it was assumed that matter only couples to one of the fields 
of the theory. The possibility of coupling matter to both fields was explored in Refs. [35- 
45]. The inclusion of other matter couplings does not seem to alleviate the problem of the 
instabilities [40, 41]. 

In this note, we address the question of whether scalar gradient instabilities auto¬ 
matically disqualifies bimetric theory or only make the theory extremely complicated phe¬ 
nomenologically since linear perturbation theory is not valid. Ref. [46] showed that the 
instabilities can be removed by taking the general relativity limit of the Hassan-Rosen the¬ 
ory, in which the effective Compton wavelength of the massive graviton goes to zero. On 
qualitative and semi-quantitative grounds, we argue for another alternative, namely that 
under the assumption of a Vainshtein mechanism, scalar instabilities have a negligible im¬ 
pact on structure formation, at least at high redshifts, and that linear structure formation 
theory can be used in the same regimes as in GR. 

We will use units where G = c = 1, where c here refers to the velocity of light. In 
the remainder of the paper, however, it will denote the relation between the asymptotic 
background metrics, as described below. 

2 The Vainshtein mechanism 

When expanding the metric fields around a flat background, higher order perturbative 
terms will become important inside the so-called Vainshtein radius, ry [47]. This regime 
does not necessarily correspond to the case where all the metric perturbations, or even 
density fluctuations, become of order unity. The existence of the Vainshtein mechanism 
for static and spherically symmetric (SSS) systems in the context of bigravity was shown 
in Refs. [48, 49]. An explicit expression of the second order term and the Vainshtein radius 
that they implied in terms of the parameters of the Hassan-Rosen theory was presented in 
Ref. [50] (see also Ref. [51]). 

Due to the complexity of the equations, the Vainshtein mechanism and possible re¬ 
covery of GR at small distances from the source have not been studied analytically for 
non-spherically symmetric and non-static systems in bimetric theory. However, in the 
following, we will assume that also for time varying density fluctuations, the Vainshtein 
mechanism will become effective at a distance similar to the corresponding Vainshtein ra¬ 
dius for the SSS-case and that GR is recovered at r <C ry. Justification of this assumption 
can be found for the DGP and cubic Galileon model in Refs. [52-54], see also Refs. [55, 56] 
for tests of the quasi-static approximation in other modified gravity models. However, for 
bimetric theory, this has not been shown explicitly and should only be regarded as an 
assumption at this point. 

We can identify the Vainshtein radius for the SSS case using the following diagonal 
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forms for and f^ v (following the notation of Ref. [48]: 

ds 2 = -F 2 dt 2 + B~ 2 dr 2 + R 2 dQ 2 , (2.1) 

ds) = -p 2 dt 2 + b 2 dr 2 + U 2 dn 2 . (2.2) 

We perturb the metrics around flat space ( g and / /t „ = c 2 g^ v ) where, in order to 
have flat backgrounds, we put 


/?o + 3/3i c + 3/3 2 c 2 + /?3 c 3 — 0, (2-3) 

-b 3/3 2 + 3/?3c + /?4C 2 = 0. (2-4:) 

c 

Note that if we instead of zero set the right hand sides of the above equations to ^FL^/m 2 , 
we obtain de Sitter solutions on the background level. Equating first and second order 
perturbation terms, we identify the Vainshtein radius as 

A/f 

(ry) 3 = qx{Pi,c) —, (2.5) 

m- 

for the metric fields X = {T 1 , B, R,p,b,U}. Here M is the mass of the source as measured 
by an observer at infinity and m the mass scale of the theory (not equal to the effective 
graviton mass). The function qx depends on the parameters of the theory and is generally 
of order unity when /% ~ c ~ 1. As an example, it has the following form for the p-field: 


c (/3i + 3 cfc + 2 c 2 f3 3 ) 

q v = --- 7T. (2.6) 

6 {j3i + 2 c/3 2 + c 2 / 3 3 ) 2 

In the following we will denote with ry the largest of the Vainshtein radii (thus setting the 
scale at which higher order terms become important), and write 


where A = m 1 and 


3 Vf o 

ry = — 2 q = MX 2 q, 
m z 


q = max.[q x (Pi,c)]. 


(2.7) 

( 2 . 8 ) 


As will be discussed in Sec. 4.4, the value of q will depend not only on the /%’s and c, but 
also on the specific metric ansatz used for g^ v and f )W . 

Numerically, we can write the Vainshtein radius as 


ry 


0.08 


/ <? 

1/3 

\ M 

f A Vl 

U 2 


Mq 

\ r Hj _ 


1/3 


kpc, 


(2.9) 


where r# = H^ 1 is the Hubble radius and h = i7o/(100km/s/Mpc). Setting A = t'h 
and M = 10 12 Mq, corresponding to a rather massive galaxy, we get ry ~ (g/h 2 ) 1 / 3 Mpc. 
Increasing the mass by a factor of a thousand (to obtain a heavy cluster), we increase the 
Vainshtein radius by a factor of ten. We thus get a radius that naturally adapts to the 
size of the object at hand. This will be of interest when studying the dynamics of galaxies 
and clusters of galaxies [51, 57]. We emphasize again that the Vainshtein radius does not 
correspond to the radius where all metric perturbations or density fluctuations become 
large. In fact, the gravitational potential <h = F 2 — 1 ~ M/r is of order <h (M/A 2 ) 2 / 3 = 
(M/r//) 2 / 3 at r = ry, which for the Sun gives < h(ry) ~ 10 -1 ' 5 and for a galaxy with mass 
M = 10 12 M q , $(ry) - 10- 7 . 
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3 Densities 


The average density within the Vainshtein radius of a spherically symmetric matter collec¬ 
tion of mass M(ry) is given by 


^ (47rr^/3) AnqX 2 ’ 

which is independent of the mass of the object. We note that this only holds for Vainshtein 
radii with a cubic form as in Eq. 2.7. For a quintic Vainshtein radius, ry oc M A 4 , the 
average density within ry decreases with mass, i.e. we will recover GR to high precision 
for small mass fluctuations and small scales but get increasing deviations for large mass 
fluctuations and large scales. For Vainshtein radii of lower order than cubic, we get large 
deviations from GR for small mass fluctuations, for which we know observationally that 
any deviation from GR has to be very small. 

In the case of A = rn = Hq 1 , 


Pv 


ZH 2 o _ 2 o 

4t xq q PcTit 


(3.2) 


where 

Pcrit = ^ ~ 2.8 • 10 11 h 2 M Q Mpc- 3 , (3.3) 

is the critical density of the Universe today. We thus expect the Vainshtein radius to 
become important for any density fluctuation larger than the critical density; a very low 
density indeed. 

Since what is relevant for the potential fluctuations are the relative density fluctuations 
5p m , defined through 

5p m Pm Pm 

Pm Pm 



we see that in order to be within the regime where the Vainshtein mechanism is active 
today, that is for 5p m > pv = (2/q)p° criv we need 6 m > 2/{VL m q), where = Pm/Pcnt- 
However, since the average density increases with decreasing scale factor, at scale factor 
a, p m = H m p[| rit a 3 , the density contrasts that we expect to be protected by the Vainshtein 
mechanism are 


<M°0 Z M°) = 


PV_ 

Pm 


2 o 3 

q£l m 


(3.5) 


4 Structure formation 

The standard picture in GR of structure formation is the following: We start out with some 
initial spectrum of perturbations. These will grow according to standard GR theory where 
the growth depends on whether the perturbations are sub- or super-horizon, whether they 
have pressure (stabilizing them on small scales) and on the background expansion. In the 
following, we only consider scalar fluctuations in pressure less matter. It is convenient to 


- 4 - 









express density fluctuations, 5 m , in terms of the amplitudes of the Fourier components, 5k, 
defined by 

5k = y j 5 m (x)exp(-ik ■ x)d 3 x, (4.1) 

since they evolve independently of each other, in the linear regime. Similarly, we denote 
the amplitudes of the Fourier components for the gravitational potential 4> by Defining 
P(k ) = (|<5fc| 2 ), where the average is over all directions, inflation generally generates initial 
conditions Pi{k ) oc k n , where n ~ 1. Since 5k oc k 2 <&k, this corresponds to Pi(&k) °c AT -4 ~ 
k~ 3 . The potential will stay roughly constant outside the horizon, preserving the shape 
of P(4>fc). When a Fourier mode re-enters the horizon at k = Ha, 

2 k 2 2 

5k ~ ~ ~3^ fc = const - (4.2) 

That is, as perturbations enter the horizon, they will be described by Ph{ k) oc fe Tl_4 ~ A: -3 . 1 
For a mode that enters the horizon during matter domination, this will take place at 
k m = H(a m )a m oc a m ' . After the entry, the mode will grow proportional to the scale 
factor, 5k oc a and the shape of the power spectrum will be modified 

„ 2 

P(k) = Pi(fc)- 5 - oc k~ 3 • k 4 = k. (4.3) 

Modes that enter earlier will only grow logarithmically, as long as radiation dominates. We 
thus expect the shape of the power spectrum for these mode to be given by 

P[k) = Pi(k) log (a rm /a r ) oc AT 3 (log k) 2 ~ k~ 3 , (4.4) 

where a rm denotes the scale factor at radiation-matter equality. Given P(k), we can 
calculate the expected variance on a given physical scale, R, as 

1 f°° 

<y\ = (\5 m (R)\ 2 ) = 7^2 J P(k)exp(-k 2 R 2 )k 2 dk, (4.5) 

where we have used a Gaussian filter to single out the Fourier modes contributing to the 
physical fluctuations. For a given scale R, a physical fluctuation 5 m (R) receives contribu¬ 
tions from all 5k where k < l/R. Note also that observationally, 5 m < 1 for R > 50 Mpc 
whereas 5k > 1 for all k ~~ 1 > 0.5 Mpc. 


4.1 Instabilities 

Imagine now that the theory of gravity at hand exhibits a gradient instability which causes 
very fast, typically exponential, growth of some (Fourier) modes, e.g. up to a scale factor 
dj, if the perturbation is within the cosmological horizon. Given an initial spectrum of 
perturbations, from Eq. 3.5 all (physical) fluctuations 5 m > 5y ~ 2a 3 /(p. m q) are protected 
from gradient instabilities by the Vainshtein mechanism, but all (sub-horizon) fluctuations 
with 5 m < 5v{di) will quickly grow to 5 m ~ 5y. Once 5 m ~ 5y, the evolution of the 

1 The power spectrum Pir(fc) is defined as the power at each k at the time when the corresponding 

Fourier mode re-enters the horizon at k = Ha. That is, it is defined at different epochs for different scales. 
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fluctuation would be equal to that of a fluctuation in GR, i.e. 5 m ~ constant during 
radiation domination and 5 m oc a during matter domination. However, with the instability 
present, the perturbation will constantly adjust to the value of the Vainshtein density, and 

2a 3 

<M a i) ~ ~ 7f^-- (4.6) 

That is, fluctuations will grow like 5 m oc af, naively in large contrast to the observed growth 
of structure. However, this predicted growth only holds for fluctuations that initially (when 
entering inside the horizon) have 5 m < 5y, and Eq. 4.6 thus only gives a lower limit on 
the amplitude of fluctuations as long as the instability is present. We next investigate 
the observational consequences of such an increased growth of the smallest fluctuation 
amplitudes, as a function of the scale factor a^ up to which the instability is present. 


4.2 Radiation domination 

At the time of photon-baryon decoupling at z ~ 1100, we have 5y ~ 5 • 10“ 9 /g and at 
matter-radiation equality at z m 3600, we have 5y ~ 1 • 10 ~ 10 /q. Let k rm correspond to 
the A:-value that crosses the horizon at a rm . 2 Since 


Hr 


Hn 


y/2Q n 

3 / 2 ” 

CLrm 


Ha 


H rrn a r 

Hn 


2fL 


V0.6 • 3600 « 46, 


(4.7) 


i.e. k rm ~ 0.015 (Mpc/li) -1 . That is, the instability will upscale all fluctuations with 
k > k rm to have a minimal value of 5 m ~ 1 • 10 ~ 10 /q at radiation-matter equality (that will 
grow to 5 m ~ 5 • 10 _7 /g today) whereas fluctuations with larger amplitude and/or size are 
not affected by the instabilities. Thus, the effect of an instability on sub-horizon matter 
perturbations during radiation domination is that for all scales smaller than approximately 
100 Mpc, we will have a minimum amplitude of fluctuations of 5 m ~ 5 • 10 /q. Using 
the Planck normalization of the amplitude of fluctuations at scale R = 8 Mpc /h [58], 
and Eq. 4.5, the average fluctuations at scales R < 100 Mpc have 5 m > 0.3. Thus, we 
can already anticipate that instabilities during radiation domination will have a negligible 
impact on the power spectrum. 


4.3 Matter domination 


If the gradient instability is present up to a scale factor a* > a rm , it will set a minimum 
value of sub-horizon perturbations at a* of 5 m (ai) > 2 a^/(Q m q). Since during matter 
domination, in GR 6 m grows proportional to the scale factor, this minimal amplitude will 
subsequently grow to a value today of 


a 0 2a/ 
mQ 


2a? 

— c m ; n ^ajj, 


n 




(4.8) 


where (5 m in(oi) is defined as the minimum perturbation value when an instability is present 
up to a scale factor a* > a rm . Here, and in the following, we have neglected the decrease 
in the growth rate when dark energy becomes dynamically important. 

2 We define horizon crossing as the point where k = Ha. 
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When constraining the maximum scale factor at which instabilities can be present, we 
need to take into account the fact that the observed 5 m depends on the observed scale R, 
and that observations corresponding to different R originate at different scale factors, a. 
During matter domination, where H ~ Hoy/tl m a ~ 3 / 2 , for sub-horizon modes 


k ^ 

#o ~ V ~a" 


(4.9) 


As noted previously, for a given scale R, a physical fluctuation 5 m (R) receives contributions 
from Fourier modes 5k where k < 1/R. If a sub-horizon instability is present for all a < a,i, 
it will affect all modes with a physical scale today of 


R 



5500 ^ 


Mpc 

~hT 


(4.10) 


Assuming sub-horizon instabilities for a < a*, we make the conservative assumption that all 
5 m (R) with R < ( Hidi ) _1 will have a minimal value today of <5 m in(oi) since the fluctuation 
will receive contributions from at least one unstable Fourier mode 5k- We can thus rule 
out all values of the scale factor a,, for which 5 m (R) < h m in(a?;), where k = R~ 3 represents 
the Fourier mode with the smallest wavelength that contributes to a fluctuation at scale 
R. This constrains the observationally viable value of a* to 


^ 6 m (R) q n m y 2 


(4.11) 


Since the observed 5 m (R) is smallest at large scales, we expect large scales to give the most 
stringent constraints on a,. However, if the Fourier modes contributing to the fluctuation 
are outside the horizon, they will not be affected by the instability. From Eq. 4.9, this 
corresponds to 


2•10~ 4 

f Mpc / 1 

k 

l h ) 


(4.12) 


We thus have two inequalities for a,i(R = k 1 ), Eqns. 4.11 and 4.12, for which if any is 
fulfilled, instabilities will not cause disagreement with current observations. 


4.4 The /3i-model 

Since we have assumed that the background geometry is flat, it is not possible to self- 
consistently derive limits with parameter values (in this case, //) corresponding to valid 
background solutions. However, regarding the asymptotic flat potential assumption as just 
a way of saying that a fluctuation should have no effect on the asymptotic space time be¬ 
haviour, we can study our solutions using parameter values that give non-flat backgrounds. 
In the case of only /3 1 , and m ~ Hq, in the de Sitter limit this corresponds to (5\ = c^ 1 = \/3 
and q = 1/3. However, this value depends on the metric ansatz Eq. 2.1. Using instead 
the metric ansatz of Ref. [50], the corresponding value is q = 28/3, that is, the Vainshtein 
radius differs by a factor of ~ 3. In the following, we will therefore account for how the 
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Figure 1. A comparison between the GR prediction for 5 m (R) (solid curved line) and the minimum 
value (5 m i n (aj) today for an instability present at a < a,i (horizontal lines) for the /3i-model using 
q = 1/3. The range of k = if” 1 for each line (indicated by the diamond) corresponds to scales that 
are sub-horizon for any a < a*, and thus are affected by the instability. 

derived limits scale with q, and keep in mind the potential uncertainty of the limits due to 
possible changes in the value of q. 

In Fig. 1, we plot the GR prediction for 5 m (R) as outlined in Sec. 4, for the /3i-model 
using q = 1/3. Note that a physical fluctuation 5 m (R) will receive contributions from 
all where k < 1/R (Eq. 4.5) and that 5 m < 1 while most 5k > 1- The horizontal 
lines correspond to the minimum perturbation amplitude today 5 m i n (a,i), for an instability 
present at a < a,i for the /3i-model. The range of k = if” 1 for each line (indicated by the 
diamond) corresponds to scales that are sub-horizon for any a < m. It is evident that the 
larger the dj, the more dramatic the modifications to the observed structures since it will 
impose a larger minimum value of 5 m over a larger range of physical scales. 

In Fig. 2, constraints given by Eqns. 4.11 and 4.12 for the /3i~model with q = 1/3 
are shown. Hatched regions indicate values not excluded by observations, since either the 
observed 5 m is larger than <5 m ; n , or all Fourier modes contributing to the fluctuation were 
outside the horizon when the instability was active. Note that this is not a constraint 
plot for dj and if” 1 , but rather for as a function of if” 1 . We can read off a minimum 
dj ss 0.015, corresponding to Zi ~ 70 from observations at if” 1 ~ 0.0015 (Mpc/h)” 1 , or 
if ~ 0.7 Gpc /h, indicated by the diamond symbol. 

However, if” 1 ~ 0.0015 (Mpc/h)” 1 corresponds to k/Ho ~ 4.5. These modes are 
only probed by CMB observations originating at z ~ 1100, thus only probing instabilities 
at d < 10” 3 , unless secondary CMB constraints such as gravitational lensing and the 
integrated Sachs-Wolfe (ISW) effect are important. Primary CMB observations do not 
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Figure 2. Constraints in the /?i-modcl with q = 1/3 on the maximum value of cq up to which 
sub-horizon gradient instabilities are effective, as a function of the observational scale R. The 
hatched regions show regions not excluded by observations, since either the observed 5 m is larger 
than <5 m i n (Eq. 4.11), or all Fourier modes contributing to the fluctuation were outside the horizon 
when the instability was active (Eq. 4.12). If any of these constraints is fulfilled, the model is in 
agreement with current observations. The white (non-hatched) region thus represents values of at in 
principle excluded. The diamond indicates the minimum allowed value of « 0.015, corresponding 
to Zi « 70 from observations at R ~ 0.7 Gpc /h. However, since low redshift observations currently 
only probe R < 100Mpc//i, the current constraint indicated by the asterisk is a, < 0.1, or Zi > 9. 


constrain instabilities at a > 10~ 3 . Low redshift galaxy observations on the other hand 
probe modes with R^ 1 > 0.01 (Mpc//i) -1 and constrain a{ < 0.1, or z\ > 9 for the /3i-model 
with q = 1/3, indicated by the asterisk in Fig. 2. Since a* oc q 1//2 (Eq. 4.11), changing the 
value of q, will change the constraint on a*. Low redshift observations thus rule out sub¬ 
horizon gradient instabilities at a scale factor larger than oq ~ 0.15Q 1 / 2 , or redshift smaller 
than Zi ~ 6.7/g 1 / 2 — 1. For example, using instead the value q = 28/3 as derived using the 
metric ansatz of Ref. [50], the corresponding limit for the /3i-model is oq < 0.5, or Zi > 1. 

We note also that as c —>• 0 or c —>• oo, it is possible that q —>■ 0 or q —>■ oo, depending 
on the values of the j 3,. As q —>• oo, the instabilities will be pushed into the far future, since 
pv ->• 0. 

5 Summary and discussion 

In this paper, we have investigated the observational consequences of sub-horizon linear 
gradient instabilities under the assumption that the Vainshtein mechanism restores GR at 
scales when the linear approximation breaks down. For static, spherically symmetric matter 
distributions, the linear approximation breaks down for density fluctuations p m > pv, 
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where pv ~ p[! ri1: for m ~ Hq. If the Vainshtein mechanism operates at similar scales also 
for quasi-static, non-spherical matter fluctuations, we have ruled out sub-horizon gradient 
instabilities at scale factors a* > 0.15(? 1 / 2 , or in terms of redshift, Zi < 6.7 /q 1 ^ 2 — 1. Here, 
q = q(Pi, c) that for , c] « 1 is of order unity, but in principle can take any value between 
zero and infinity. For the /3i-model with all /% = 0 except [3\ , depending on the assumed 
value of q, we can constrain > 1 — 10. This is to be compared to theoretical predictions 
for the same model, showing that the instability will persist down to a redshift of Zi ~ 0.5 
[22], At face value, this indicates a serious tension between observational constraints and 
theoretical predictions for the /3i-model. However, the difference is sufficiently small that 
the final verdict on the (3 \-model needs to await a more detailed analysis. For cases where 
the theoretical prediction and the observational constraints differ by larger amounts, the 
method outlined in this paper should be sufficient to rule in or out the models, in terms of 
the effects of gradient instabilities. 

Since the largest observable scales are the most efficient in constraining the effect 
of instabilities, we note that a study of secondary effects on the CMB, particularly the 
integrated Sachs-Wolfe (ISW) effect, could potentially yield more stringent constraints on 
models. The ISW effect probes time variations in the gravitational potentials between 
the CMB decoupling and the present epoch. In a matter dominated universe, we expect 
potentials to be constant and the ISW effect to be zero. If the effect from curvature and/or 
dark energy on the background expansion is non-negligible, potentials will decay with time 
and we will observe slightly higher CMB temperatures in the direction of matter over 
densities. This effect has been detected using the angular cross-power spectrum between 
the CMB temperature anisotropies and large-scale structure at around 3cr significance, see 
[28] for a list of references and a study of the observational effects on a bimetric model not 
affected by gradient instabilities. If density fluctuations would grow exponentially, so would 
the potentials. The potential will only vary in time exponentially during a very short epoch 
directly after entering inside the horizon, after which it will grow as oc a 2 . Since what is 
observed is decaying potentials, ISW observation can potentially put stringent constraints 
on models plagued by gradient instabilities. It also remains to be investigated to what 
extent the qualitative analysis presented here can also be applied to the tensor sector. 

Acknowledgments 

We would like to thank Yashar Akrami, Adam R. Solomon, Frank Konnig, Fawad Hassan 
and Angnis Schmidt-May for many interesting discussions in relation to this work. E.M. 
acknowledges support for this study by the Swedish Research Council. 

References 

[1] M. Fierz and W. Pauli, On relativistic wave equations for particles of arbitrary spin in an 
electromagnetic field , Proc.Roy.Soc.Lond. A173 (1939) 211 -232. 

[2] C. Isham, A. Salam, and J. Strathdee, F-dominance of gravity, Phys.Rev. D3 (1971) 867-873. 


- 10 - 


[3] D. Boulware and S. Deser, Can gravitation have a finite range?, Phys.Rev. D6 (1972) 
3368-3382. 

[4] S. F. Hassan and R. A. Rosen, Bimetric gravity from ghost-free massive gravity, Journal of 
High Energy Physics 2 (Feb., 2012) 126, [arXiv : 1109.3515]. 

[5] S. F. Hassan and R. A. Rosen, Resolving the Ghost Problem in Nonlinear Massive Gravity , 
Physical Review Letters 108 (Jan., 2012) 041101, [arXiv: 1106.3344]. 

[6] C. de Rham and G. Gabadadze, Generalization of the Fierz-Pauli action, Physical Review D 
82 (Aug., 2010) 044020, [arXiv: 1007.0443], 

[7] C. de Rham, G. Gabadadze, and A. J. Tolley, Resummation of Massive Gravity, Physical 
Review Letters 106 (June, 2011) 231101, [arXiv: 1011.1232]. 

[8] S. Hassan and R. A. Rosen, Resolving the Ghost Problem in non-Linear Massive Gravity, 
Phys.Rev.Lett. 108 (2012) 041101, [arXiv: 1106.3344], 

[9] S. Hassan, R. A. Rosen, and A. Schmidt-May, Ghost-free Massive Gravity with a General 
Reference Metric, JHEP 1202 (2012) 026, [arXiv: 1109.3230]. 

[10] S. Hassan, A. Schmidt-May, and M. von Strauss, Proof of Consistency of Nonlinear Massive 
Gravity in the Stuckelberg Formulation, Phys.Lett. B715 (2012) 335-339, 

[arXiv: 1203.5283], 

[11] M. von Strauss, A. Schmidt-May, J. Enander, E. Mortsell, and S. Hassan, Cosmological 
Solutions in Bimetric Gravity and their Observational Tests, JCAP 1203 (2012) 042, 

[arXiv: 1111.1655], 

[12] M. S. Volkov, Cosmological solutions with massive gravitons in the bigravity theory, JHEP 
1201 (2012) 035, [arXiv: 1110.6153], 

[13] D. Comelli, M. Crisostomi, F. Nesti, and L. Pilo, FRW Cosmology in Ghost Free Massive 
Gravity, JHEP 1203 (2012) 067, [arXiv: 1111.1983], 

[14] Y. Akrami, T. S. Koivisto, and M. Sandstad, Accelerated expansion from ghost-free bigravity: 
a statistical analysis with improved generality, Journal of High Energy Physics 3 (Mar., 2013) 
99, [arXiv: 1209.0457], 

[15] Y. Akrami, T. S. Koivisto, and M. Sandstad, Cosmological constraints on ghost-free 
bigravity: background dynamics and late-time acceleration, arXiv : 1302.5268. 

[16] H. Nersisyan, Y. Akrami, and L. Amendola, Consistent metric combinations in cosmology of 
massive bigravity, arXiv: 1502.0398. 

[17] D. Comelli, M. Crisostomi, and L. Pilo, Perturbations in Massive Gravity Cosmology, JHEP 
1206 (2012) 085, [arXiv: 1202.1986], 

[18] N. Khosravi, H. R. Sepangi, and S. Shahidi, Massive cosmological scalar perturbations, 
Phys.Rev. D86 (2012) 043517, [arXiv: 1202.2767], 

[19] M. Berg, I. Buchberger, J. Enander, E. Mortsell, and S. Sjors, Growth Histories in Bimetric 
Massive Gravity, JCAP 1212 (2012) 021, [arXiv : 1206.3496]. 

[20] M. Fasicllo and A. J. Tolley, Cosmological Stability Bound in Massive Gravity and Bigravity, 
JCAP 1312 (2013) 002, [arXiv: 1308.1647], 

[21] F. Konnig, A. Patil, and L. Amendola, Viable cosmological solutions in massive bimetric 
gravity, JCAP 1403 (2014) 029, [arXiv: 1312.3208], 


- 11 - 


[22] F. Konnig and L. Amendola, Instability in a minimal bimetric gravity model, Phys.Rev. D90 
(2014) 044030, [arXiv: 1402.1988], 

[23] D. Comclli, M. Crisostomi, and L. Pilo, FRW Cosmological Perturbations in Massive 
Bigravity , Phys.Rev. D90 (2014) 084003, [arXiv: 1403.5679]. 

[24] A. De Felice, A. E. Gumrukgiioglu, S. Mukohyama, N. Tanahashi, and T. Tanaka, Viable 
cosmology in bimetric theory , JCAP 1406 (2014) 037, [arXiv : 1404.0008]. 

[25] A. R. Solomon, Y. Akrami, and T. S. Koivisto, Linear growth of structure in massive 
bigravity , JCAP 1410 (2014) 066, [arXiv: 1404.4061]. 

[26] F. Konnig, Y. Akrami, L. Amendola, M. Motta, and A. R. Solomon, Stable and unstable 
cosmological models in bimetric massive gravity , Phys.Rev. D90 (2014) 124014, 

[arXiv: 1407.4331], 

[27] M. Lagos and P. G. Ferreira, Cosmological perturbations in massive bigravity, JCAP 1412 
(2014) 026, [arXiv: 1410.0207], 

[28] J. Enander, Y. Akrami, E. Mortsell, M. Renneby, and A. R. Solomon, Integrated Sachs-Wolfe 
effect in massive bigravity , Phys.Rev. D91 (2015) 084046, [arXiv : 1501.0214]. 

[29] F. Konnig, Higuchi Ghosts and Gradient Instabilities in Bimetric Gravity , Phys.Rev. D91 
(2015) 104019, [arXiv: 1503.0743], 

[30] K. Aoki, K.-i. Maeda, and R. Namba, Stability of the Early Universe in Bigravity Theory, 
arXiv:1506.0454. 

[31] G. Cusin, R. Durrer, P. Guarato, and M. Motta, Gravitational waves in bigravity cosmology, 
JCAP 1505 (2015), no. 05 030, [arXiv: 1412.5979], 

[32] M. Johnson and A. Terrana, Tensor Modes in Bigravity: Primordial to Present, 
arXiv:1503.0556. 

[33] L. Amendola, F. Konnig, M. Martinclli, V. Pettorino, and M. Zumalacarregui, Surfing 
gravitational waves: can bigravity survive growing tensor modes?, JCAP 1505 (2015), no. 05 
052, [arXiv: 1503.0249], 

[34] M. Fasicllo and R. H. Ribeiro, Mild bounds on bigravity from primordial gravitational waves, 
arXiv:1505.0040. 

[35] Y. Akrami, T. S. Koivisto, and A. R. Solomon, The nature of spacetime in bigravity: Two 
metrics or none?, General Relativity and Gravitation 47 (Jan., 2015) 1838, 

[arXiv: 1404.0006], 

[36] C. de Rham, L. Heisenberg, and R. H. Ribeiro, On couplings to matter in massive 
(bi-)gravity, Class. Quant. Grav. 32 (2015) 035022, [arXiv : 1408.1678]. 

[37] C. de Rham, L. Heisenberg, and R. H. Ribeiro, Ghosts and matter couplings in massive 
gravity, bigravity and multigravity, Phys.Rev. D90 (2014) 124042, [arXiv : 1409.3834], 

[38] J. Enander, A. R. Solomon, Y. Akrami, and E. Mortsell, Cosmic expansion histories in 
massive bigravity with symmetric matter coupling , JCAP 1501 (2015) 006, 

[arXiv: 1409.2860], 

[39] A. Schmidt-May, Mass eigenstates in bimetric theory with matter coupling, JCAP 1501 
(2015) 039, [arXiv: 1409.3146], 

[40] D. Comclli, M. Crisostomi, K. Koyama, L. Pilo, and G. Tasinato, Cosmology of bigravity 
with doubly coupled matter , JCAP 1504 (2015) 026, [arXiv : 1501.0086]. 


- 12 - 


[41] A. E. Giimriikguoglu, L. Heisenberg, S. Mukohyama, and N. Tanahashi, Cosmology in 
bimetric theory with an effective composite coupling to matter, JCAP 1504 (2015), no. 04 
008, [arXiv: 1501.0279], 

[42] Y. Akrami, T. S. Koivisto, D. F. Mota, and M. Sandstad, Bimetric gravity doubly coupled to 
matter: theory and cosmological implications, JCAP 1310 (2013) 046, [arXiv : 1306.0004]. 

[43] A. R. Solomon, J. Enander, Y. Akrami, T. S. Koivisto, F. Knnig, et al., Cosmological 
viability of massive gravity with generalized matter coupling, JCAP 1504 (2015), no. 04 027, 
[arXiv: 1409.8300], 

[44] A. Emir Gmrkolu, L. Heisenberg, and S. Mukohyama, Cosmological perturbations in massive 
gravity with doubly coupled matter, JCAP 1502 (2015) 022, [arXiv : 1409.7260]. 

[45] L. Heisenberg, More on effective composite metrics, arXiv : 1505 . 0296. 

[46] Y. Akrami, S. Hassan, F. Konnig, A. Schmidt-May, and A. R. Solomon, Bimetric gravity is 
cosmologically viable, arXiv: 1503.0752. 

[47] A. I. Vainshtein, To the problem of nonvanishing gravitation mass, Physics Letters B 39 
(May, 1972) 393-394. 

[48] M. S. Volkov, Hairy black holes in the ghost-free bigravity theory, Phys.Rev. D85 (2012) 
124043, [arXiv: 1202.6682], 

[49] E. Babichev and M. Crisostomi, Restoring general relativity in massive bigravity theory, 
Phys.Rev. D88 (2013), no. 8 084002, [arXiv: 1307.3640], 

[50] J. Enander and E. Mortsell, Strong lensing constraints on bimetric massive gravity, Journal 
of High Energy Physics 10 (Oct., 2013) 31, [arXiv : 1306.1086]. 

[51] J. Enander and E. Mortsell, On stars, galaxies and black holes in massive bigravity, 

arXiv:1507.0091. 

[52] F. Schmidt, Self-consistent cosmological simulations of DGP braneworld gravity, Physical 
Review D 80 (Aug., 2009) 043001, [arXiv:0905.0858]. 

[53] H. A. Winther and P. G. Ferreira, The Vainshtein mechanism beyond the quasi-static 
approximation, arXiv: 1505.0353. 

[54] R. Brito, A. Terrana, M. Johnson, and V. Cardoso, Nonlinear dynamical stability of infrared 
modifications of gravity, Phys.Rev. D90 (2014) 124035, [arXiv : 1409.0886]. 

[55] C. Llinares and D. F. Mota, Releasing Scalar Fields: Cosmological Simulations of 
Scalar-Tensor Theories for Gravity Beyond the Static Approximation, Physical Review 
Letters 110 (Apr., 2013) 161101, [arXiv: 1302.1774], 

[56] S. Bose, W. A. Hellwing, and B. Li, Testing the quasi-static approximation in f(R) gravity 
simulations, Journal of Cosmology and Astroparticle Physics 2 (Feb., 2015) 34, 

[arXiv: 1411.6128], 

[57] K. Koyama and J. Sakstein, Astrophysical Probes of the Vainshtein Mechanism: Stars and 
Galaxies, arXiv: 1502.0687. 

[58] Planck Collaboration, P. A. R. Ade, N. Aghanim, M. Arnaud, M. Ashdown, J. Aumont, 

C. Baccigalupi, A. J. Banday, R. B. Barreiro, J. G. Bartlett, and et al., Planck 2015 results. 
XIII. Cosmological parameters, ArXiv e-prints (Feb., 2015) [arXiv : 1502.0158]. 


- 13 - 


